——— Introducción al análisis espacial ———
Serie de talleres en “Herramientas de análisis cuantitativo y su
aplicación en la conservación de la biodiversidad”
Serie de talleres en “Herramientas de análisis cuantitativo y su
aplicación en la conservación de la biodiversidad”
Este taller ha sido desarrollado por:
- Juan Zuloaga.
- Prof. Melissa Guzmán (University of Southern California).
- Keyvette Tabb (University of Southern California)
1. Objetivos de aprendizaje
Al final de este taller usted podrá:
- Identificar las características de los datos espaciales, en sus dos formatos: ráster y vectorial.
- Conocer los paquetes y funciones en R más utilizados para el manejo de los dos formatos.
- Utilizar bases de datos de libre acceso disponibles en la web.
- Integrar y graficar los dos formatos
- Implemetar algunos análisis espaciales que pueden servir de herramientas para la toma de decisiones en la conservación de la biodiversidad.
2. ¿Qué hacer si no tienes conocimiento previo en R?
Si no tiene concocimientos de R le sugerimos tomar estos cursos que le ayudarán a tener bases y desarrollar sus habilidades.
3. Paquetes en R que vamos a utilizar
terra: Análisis espacial de datos (Hijmans et al 2022).raster: Análisis y modelamiento de datos geográficos (Hijmans et al 2022).sf: Simple Features (en inglés) permite el manejo de datos vectoriales (Pebesma et al 2022).stars: permite el manejo de colecciones espacio-temporales y cubos de datos en formato vectorial y raster (Pebesma et al 2022).gdalcubes: permite el manejo de cubos de datos, por ejemplo observaciones de la tierra a través de colecciones de imágenes de satélite (Appel et al 2022).geodataPaquete para decargar datos geogáficos (Hijmans et al 2022).microbenchmark: Mide el tiempo de procesamiento (parece ser mejor que la funciónsystem.time()que viene en el paquete base). (microbenchmark)
raster y terra son paquetes en R que
permiten el manejo de datos en formato ráster. terra ha
sido desarrollado para reemplazar raster, es mucho más
rápido y posee mas funcionalidas (entre ellas el manejo de datos
vectoriales). En
este
vínculo puede ver las principales funciones de terra y
la comparación de las principales funciones que se encontraban en
raster y cómo se llaman ahora en terra.
Por otro lado, el paquete sf (simple features, en
inglés) ha reemplazado completamente el paquete sp. Vamos
entonces a enfocarnos en sf que permite el manejo de datos
vectoriales, como por ejemplo: puntos, líneas y polígonos.
Al final veremos algunas las aplicaciones de los paquetes
stars y gdalcubes que permiten el manejo de
cubos de datos; como por ejemplo imágenes de satélite que es una
colección de varias bandas y años.
Antes de empezar, vamos a instalar (si es necesario) y cargar estos
paquetes utilizando la función library().
4. Para empezar
4.1. Configurar folder de trabajo
Inicie RStudio y abra un nuevo R Script utilizando el menu: File > New File > R Script. Usted podrá copiar los códigos que aquí presentamos y probarlos localmente en su computador. Además habrá unos ejercicios que le permitirán consolidar y aumentar sus habilidades en desarrollar códigos y realizar análisis espaciales.
Digámos que crea un directorio llamado “analisis_espacial” para guardar los archivos de trabajo.
Recuerde que cuando abre RStudio y un nuevo R Script, es importante
saber donde esta trabajando con respecto a su disco duro. Para eso puede
utilizar la función getwd().
Para asignar el nuevo directorio de trabajo (working directory, en
inglés), utilice la función setwd().
# getwd() # consigue el directorio en el cual está trabajando
# setwd("C:/analisis_espacial")
# en mi caso voy es
# setwd("C:/Talleres_R/IntroR_espacial")Ahora que hemos asignado nuestro directorio de trabajo personal, R sabe donde buscar para cargar o guardar archivos!
4.2. Configurar folder de archivos temporales
Cuando se manipulan o se hacen análisis con rásters en R, se generan archivos temporales.
Generalmente uno no sabe esto y menos en donde se guardan estos archivos.
Utilice esta función para ver por defecto en donde se guardan:
Ahora bien, para tener más control de los archivos temnporales es buena practica configurar un folder en su proyecto que se encargue de guardar estos archivos. Al final de su proyecto los podrá borrar. Veremos su utilidad más adelante.
rasterOptions(tmpdir = "./Temp_para_borrar") #Note que estoy utilizando la ruta relativa para crear un folder en mi proyecto "C:/Talleres_R/IntroR_espacial"
rasterTmpFile()## [1] "C:/Talleres_R/IntroR_espacial/Temp_para_borrar/r_tmp_2024-07-21_194728.352313_28580_92047.grd"
Configurar otros folders para su proyecto
Esto es opcional, pero le permitirá manejar mejor su proyecto.
- Es bueno crear un folder para guarar todos los datos que bajamos de la red. Por ejemplo creemos un folder llamado “input_data”.
if(!dir.exists("C:/Talleres_R/IntroR_espacial/Input_data")){
dir.create("C:/Talleres_R/IntroR_espacial/Input_data")
}else{
print("Folder Input_data ya existe")
}## [1] "Folder Input_data ya existe"
- Tambien otro folder para guadar sus rasters intermedias que genera dentro de su proyecto pero que probablemente no son los resultados finales.
if(!dir.exists("C:/Talleres_R/IntroR_espacial/Intermediate")){
dir.create("C:/Talleres_R/IntroR_espacial/Intermediate")
}else{
print("Folder Intermediate ya existe")
}## [1] "Folder Intermediate ya existe"
- Finalmente si quiere un folder que contenga los raster finales.
if(!dir.exists("C:/Talleres_R/IntroR_espacial/Outputs")){
dir.create("C:/Talleres_R/IntroR_espacial/Outputs")
}else{
print("Folder Outputs ya existe")
}## [1] "Folder Outputs ya existe"
5. ¿Por qué es importante sabe manejar y analizar información espacial?
- Entender patrones de biodiversidad
- Hacer prediciones (como porjemplo sobre la distribución de las espcies y el efecto de cambio climático).
- Planear areas de monitoreo
- Planeación de la conservación
- Entre otras
6. ¿Cuáles son los formatos en que encontramos información espacial?
Los principales formatos en el cual encontramos información espacial son: vectorial y ráster.
Los datos vectoriales pueden representar puntos, líneas o polígonos. En general, este formato es común para almacenar y representar observaciones (ej., presencia de especies), objetos lineales (ej., como vías) y áreas geográficas (ej., límites administrativos).
Los datos en ráster (que se concocen también como malla o grilla) almacenan pixeles (celdas). Imagine una fotografía en su cámara digital, o las imágenes de su televisor, o las obtenidas por los satélites, etc. Estas imágenes están compuestas de pixeles.
En adelante utilizaremos el término técnico ráster para
hacer referencia a este formato de datos, pero para dar fluidez al texto
tambien puede encontrar malla o grilla.
Utilizaremos también las palabras pixel o
celda para hacer referencia a la unidad mínima de un ráster
que contiene información.
7. Ráster
Como mencionamos, el ráster es un formato de datos que almacena imágenes digitales en forma de pixeles (o celdas). En general es una cuadrícula con varias celdas que contienen valores.
Empecemos con unos ejemplos para visualizar y entender mejor de lo que estamos hablando.
7.1. Creando rásters
Vamos a crear un objecto ráster de 10 columnas y 10 filas, utilizando
la funcion rast() del paquete terra.
## class : SpatRaster
## dimensions : 10, 10, 1 (nrow, ncol, nlyr)
## resolution : 36, 18 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
Al explorar el objeto r1 vemos lo siguiente:
- el parámetro
class:esSpatRaster(es decir, es básicamente un ráster). - el parámetro
dimensions:es de 10 filas por 10 columnas (tal y como la definimos) y una capa (nlyr).
Pero además vemos que la función le ha asignado otros atributos como:
cood. ref.:la función asigna por defecto coordenadas geográficas lon/lat, usando WGS 84 (World Geodetic System, por sus siglas en inglés).
Ahora bien, la función ha generado únicamente la estructura del ráster, pero todavía el ráster no tiene valores.
Asignemos valores a cada pixel con la funcion values().
Note que la función ncell() calcula el número de pixeles en
el ráster. En este caso vamos asignar valores secuenciales.
## class : SpatRaster
## dimensions : 10, 10, 1 (nrow, ncol, nlyr)
## resolution : 36, 18 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## name : lyr.1
## min value : 1
## max value : 100
La función values() ha asignado varias cosas:
- los valores de 1 a 100 (
min value:ymax value:, respectivamente) extent:por defecto la extensión es en coordenadas geográficas. En este caso longitud (-180 a 180) y latitud (-90 a 90; de polo sur a polo norte).resolution:es la resolución de la grilla, es decir el tamaño de cada pixel o celda, se calcula con base en el número de filas y columnas que indicamos al principio (dimensions:) y con base en elextent:. Es decir, como definimos 10 por 10, la función asigna el siguiente tamaño de pixel 36 x 18 (grados de longitud(x) y latitud(y), respectivamente).- el nombre del ráster (
name:lyr.1). Lo podemos cambiar más adelante con un nombre que refleje lo que contiene el ráster. - el parámetro
source:indica la fuente del ráster. En este caso el objeto está en la memoria virtual del computador. Más adelante podrá ver que si descarga un ráster desde su computador podrá utilizar la funciónsources()para conocer la ruta donde está guardado el ráster.
Visualicemos el ráster que hemos creado, utlizando la función
plot() y la función text() para agregar los
valores de cada pixel sobre el ráster.
plot(r1) # uilizamos la función plot para visualizar r1
text(r1) # utilizamos la función text para adicionar los valores de cada celdaLa grilla que se genera posee valores enteros y la función
plot() genera la escala de la barra vertical.
Ya hemos creado nuestro primer ráster!
7.2. Cobertura espacial (extent)
extent: de un ráster representa la cobertura geográfica
del ráster, incluyendo celdas sin valor (NAs).
Creemos un ráster basados en el objeto r1 pero
reemplacemos los valores de la primera y segunda filas con
NA (no values).
## class : SpatRaster
## dimensions : 10, 10, 1 (nrow, ncol, nlyr)
## resolution : 36, 18 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## name : lyr.1
## min value : 21
## max value : 100
Al visualizar el objeto vemos que las primera dos filas del objeto no
se observan (tal y como la creamos), pero la extensión espacial es la
misma
extent: -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
plot(r_NA, legend=F)
text(r_NA, cex=0.9) # utilizamos la función text para adicionar los valores de cada celdaComparemos los dos objetos r1 y r_NA.
Como puede ver los dos rásters tiene en mismo extent
pero en el segundo no hay datos en las primeras dos filas
Aunque parezca trivial es muy importante tener esto en cuenta, más adelante veremos algunos efectos de esto cuando aplicamos métodos de análisis.
par(mfrow=c(1,2))
# r1
plot(r1, legend=F)
text(r1, cex=0.9)
# Adicionemos un punto y las coordenadas en cada esquina de la extensión espacial.
points(c(xmin(r1), xmax(r1), xmin(r1), xmax(r1)), c(ymax(r1), ymax(r1), ymin(r1), ymin(r1)), pch=20, col=c("blue", "red"),cex=4)
# Adicionemos texto
graphics::text(-40, 97, "xmin=-180; ymax=90", pos=2, col="blue", cex=1.5)
graphics::text(180, 97, "xmax=180; ymax=90", pos=2, col="red", cex=1.5)
graphics::text(-40, -97, "xmin=-180; ymin=-90", pos=2, col="blue", cex=1.5)
graphics::text(180, -97, "xmax=180; ymin=-90", pos=2, col="red", cex=1.5)
# r_NA
plot(r_NA, legend=F)
text(r_NA, cex=0.9)
# Adicionemos un punto y las coordenadas en cada esquina de la extensión espacial.
points(c(xmin(r_NA), xmax(r_NA), xmin(r_NA), xmax(r_NA)), c(ymax(r_NA), ymax(r_NA), ymin(r_NA), ymin(r1)), pch=20, col=c("blue", "red"),cex=4)
# Adicionemos texto
graphics::text(-40, 97, "xmin=-180; ymax=90", pos=2, col="blue", cex=1.5)
graphics::text(180, 97, "xmax=180; ymax=90", pos=2, col="red", cex=1.5)
graphics::text(-40, -97, "xmin=-180; ymin=-90", pos=2, col="blue", cex=1.5)
graphics::text(180, -97, "xmax=180; ymin=-90", pos=2, col="red", cex=1.5)7.3. Utilidad de crear un ráster ficticio
El ráster que hemos creado es lo que conocemos como un ráster
ficticio (dummy raster en inglés). Las utilidades de
construir un ráster ficticio son varias, por ejemplo:
podemos utilizar el ráster ficticio para probar funciones/operaciones y verificar los resultados antes de aplicarlas a un ráster grande y con muchos datos. Es decir, podemos verificar que nuestras funciones arrojen los resultados esperados, antes de aplicar dichas funciones a bases de datos reales.
podemos compatir en la web nuestros códigos para resolver un problema al cual no le encontramos solución. Es decir compartimos un ejemplo reproducible de lo que queremos hacer y así otras personas pueden sugerir modificaciones. Unas lecturas de la importancia de crear un buen ejemplo reproducible: Elio Campiteli, stackoverplow.com.
Me voy a explicar con un ejemplo para el primer caso.
Primero creemos otro ráster, por ejemplo: r2 será igual
a r1+1.
Revise el gráfico de abajo. ¿Es lo esperado?
r2 <- r1+1 # r2 es el resultado de sumarle 1 a r1 creada anteriormente
par(mfrow=c(1,2)) #definamos una ventana de una fila y dos columnas para presentar nuestros rásters
plot(r1, legend=F, main="r1")
text(r1)
plot(r2, legend=F, main="r2")
text(r2)Creemos ahora otro ráster r3 con base en r1
y r2.
La formula que vamos a utilizar es (r1+r2)/10. Es decir, sumar los valores de r1 y r2 y luego dividirlos entre 10.
En el gráfico de abajo se encuentran los resultados.
¿Es lo esperado?
Parece que sí. Por ejemplo, la primera celda es 0.3, que es el resultado de (1+2)/10 = 0.3
Este es un ejemplo muy simple, pero espero que muestre la utilidad de crear rásters ficticios.
7.4. Apilando rásters
Una de las funcionalidades más interesantes con rásters es que las podemos apilar, es decir crear un objeto con un conjunto de capas .
Por ejemplo, apilemos r1, r2,
r3 en un objeto que llamaremos r_stack,
utilizando la funcion c() del paquete terra
(la función análoga en el paquete raster es
stack() o brick()).
## class : SpatRaster
## dimensions : 10, 10, 3 (nrow, ncol, nlyr)
## resolution : 36, 18 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## names : lyr.1, lyr.1, lyr.1
## min values : 1, 2, 0.3
## max values : 100, 101, 20.1
Vemos que tenemos ahora un objeto con tres capas que se pueden
identificar por el nombre (name:en este caso los nombres
son iguales) con sus respectivos valores: min values y
max values:.
Dentro de este objeto podemos visualizar cada capa de varias formas:
- La pila de rásters es como una lista de capas, entonces usted puede
acceder a cada una utilizando el nombre del objeto seguido de la capa
que quiere entre dos parentesis cuadrados, por ejemplo la capa 2
r_stack[[2]].
## class : SpatRaster
## dimensions : 10, 10, 1 (nrow, ncol, nlyr)
## resolution : 36, 18 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## name : lyr.1
## min value : 2
## max value : 101
- La otra forma es utilizando el signo pesos
$despues del objeto (r_stack$) y le aparecera la lista de las capas que se encuentran en el objeto, identificadas con el nombre (name:). En este caso selecionemos la primera capa (r_stack$lyr.1). El problema es que todas las capas tienen el mismo nombre, entonces usted obtendrá este error:
## Error: [subset] you cannot select a layer with a name that is not unique
Aprovechemos y cambiemos los nombres de las capas en esta pila.
names(r_stack) <- c("capa_1", "capa_2", "capa_3") # encadenamos (c) el nombre de cada capa utilizando comillas.
r_stack## class : SpatRaster
## dimensions : 10, 10, 3 (nrow, ncol, nlyr)
## resolution : 36, 18 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## names : capa_1, capa_2, capa_3
## min values : 1, 2, 0.3
## max values : 100, 101, 20.1
Ahora sí intentemos utilizando el signo $. Parece que si
funciona.
## class : SpatRaster
## dimensions : 10, 10, 1 (nrow, ncol, nlyr)
## resolution : 36, 18 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## name : capa_2
## min value : 2
## max value : 101
7.5. Resolución espacial
La resolución espacial de un ráster es básicamente el tamaño del pixel que a su vez representa el área sobre la superficie terrestre.
Por ejemplo en las imágenes Landsat 8 de la (NASA) se muestra la ciudad de Reykjavik en Islandia (July 7, 2019) en donde se pueden apreciar tres resoluciones (30m, 100m y 300m). Es decir desde una resolución fina con más detalle (cada pixel en el terreno es de 30m por 30m) a una más gruesa y borrosa (cada pixel en el terreno es de 300m por 300m).
Como la resolución más fina tiene más detalle, es decir má información, se podrían diferenciar los elementos del terreno y hacer una clasificación detallada y precisa de la superficie terrestre. En la más gruesa es un poco más difícil y se tendría que hacer unas generalizaciones sobre lo que creemos que predomina en cada pixel.
Veamos ahora el concepto de resolución con nuestras rásters
ficticias. Creemos otro ráster r4 de 100 columnas y 100
filas y adicionemos el nombre “capa_4”.
r4 <- rast(ncol=100,
nrow=100,
)
values(r4) <- 1:ncell(r4)
names(r4) <- "capa_4" # Note que el objeto es r4 y el nombre de la capa lo cambiamos a "capa_4"Comparemos r4 con r_stack$capa_1.
## class : SpatRaster
## dimensions : 100, 100, 1 (nrow, ncol, nlyr)
## resolution : 3.6, 1.8 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## name : capa_4
## min value : 1
## max value : 10000
## class : SpatRaster
## dimensions : 10, 10, 1 (nrow, ncol, nlyr)
## resolution : 36, 18 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## name : capa_1
## min value : 1
## max value : 100
La grilla r4 es mucho mas fina que
r_stack$capa_1 porque:
la resolución de
r4es de 3.6 por 1.8. Es decir que cada pixel es un rectangulo con un lado de 3.6 y el otro 1.8.y en el ráster
r_stack$capa_1los pixeles son de mayor tamaño (rectangulares tambien) de 36 por 18.
Note que las celdas no necesariamente tienen que ser cuadradas,
especialmente cuando tenemos coordenadas geográficas. Más adelante
veremos como los rásters projectados sí tienen celdas cuadradas (bueno,
en lo posible). Es decir x, y en resolution: deberían tener
valores iguales.
Volviendo al ejemplo anterior, podemos también decir que la
resolución de r_stack$capa_1 es mucho más gruesa. de hecho
10 veces más gruesa que r4.
Como la resolución de r4 es más fina, entonces hay mas
celdas (10,000; dimensions: 100 por 100) que en
r_stack$capa_1 (100; dimensions: 10 por
10).
Note que el área seleccionada (extent:) y las
coordenadas (coord. ref.:) son las mismas. Visualicemos las
dos rásters para ver sus diferencias:
par(mfrow=c(1,2))
plot(r4, legend=F, main="r4, resolución fina")
plot(r_stack$capa_1, legend=F, main="r_stack$capa_1, resolución gruesa")Ahora bien, ¿Podríamos crear una pila de rásters con área
(extent) y coordenadas (coord. ref.)
similares, pero resoluciones y dimensiones diferentes?
Veamos:
## Error: [rast] number of rows and/or columns do not match
Pues NO, el Error nos dice que el número de columnas y
filas (dimensions) no coinciden.
7.6. ¿Qué más podemos hacer con rásters?
Ya vimos algunas cosas que podemos hacer con rásters, como por ejemplo apilarlas. Pero hay algunos métodos importantes cuando se quieren hacer análisis utilizando rásters. Cuando nos referimos a métodos hacemos referencia a una serie de funciones que se pueden aplicar a uno o varios rásters para obtener información relevante en su investigación.
No podemos en este taller ver todas las funciones, pero al final usted tendrá la capacidad de explorar y aplicar más funciones.
Por ejemplo, ya vimos cómo podemos crear un ráster utilizando la
función rast(), apilar varios rásters con c(),
y cómo seleccionar algunas capas dentro de una pila de rásters
utilizando [[]] o $.
Ahora veamos algunos métodos y funciones que puden ser mucha utilidad en el análisis espacial.
7.6.1. Métodos locales
Los métodos locales hacen referencia a hacer cálculos entre celdas que se sobreponen en una pila de rásters. Es decir la localización espacial de las celdas es la misma. La figura muestra cómo podríamos calcular la media de dos rásters que se sobreponen.
Hay varias funciones para hacer cálculos de este tipo. Por ahora
utilizaremos la función app() para demostrar este
método.
Veamos el código en R para calcular la media de la pila de tres capas
que hemos creado; es decir la funcion app() en
terra toma los tres valores de cada pixel de cada capa que
se sobreponen y calcula la media (mean). Esto se hace en
toda la extensión del área que definimos.
r_stack_media <- app(x = r_stack, # Este es el objeto que representa la pila de rásters
fun = mean # Esta es la función para calcular la media
)
r_stack_media## class : SpatRaster
## dimensions : 10, 10, 1 (nrow, ncol, nlyr)
## resolution : 36, 18 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## name : mean
## min value : 1.1
## max value : 73.7
Note que el nuevo objeto r_stack_media representa la
media de las tres capas en r_stack. R asigna el nombre de
la función a la capa final en esta caso name: es
mean.
Visualicemos el promedio de esto tres rásters.
Tips!
Note que los decimales están indicados con un punto .,
para los miles se utiliza la coma ,. Los países del norte
utilizan esta nomenclatura y R tambien.
En Colombia es al contrario (no sabría decir si en toda latinoamérica es igual). Por qué sucede esto? No lo sé.
Lo que sé es que usted tendrá que ajustar las bases de datos para que tengan el mismo formato requerido en R. En especial tablas en excel que le gustaría incorporarlas a un análisis espacial.
Usted ya puede ver las potenciales aplicaciones de estas funciones en un pila de rásters. Por ejemplo, calcular la media o la varianza de la temperatura durante un año (una pila de 12 capas) en su área de estudio.
Ahora veamos un ejemplo utilizando datos globales reales!!!!
Datos globales y algunas aplicaciones
Vamos ahora a manejar algunos bases de datos de acceso gratuito (open access, en inglés) que pueden ser utilizadas en análisis espaciales.
Por ejemplo, los datos globales de clima, como temperatura y precipitación, pueden ser utilizados ampliamente para entender los efectos del cambio climático en la distribución de las especies.
WorldClim y CHELSA son dos iniciativas que ofrecen datos de buena calidad a nivel global.
Existen además algunos paquetes en R que permiten descargar esta
bases de datos. Por ejemplo, en el paquete raster existe la
función getData() que permite descargar datos de WorldClim.
Sin embargo, esta función será reemplazada por el paquete
geoData. Veamos un ejemplo
Temperatura media anual a nivel global (WorldClim)
Descarguemos la temperatura media anual a nivel global a una resolución gruesa (10 minutos que es ~340 km2) desde WorldClim. Tenga en cuenta que descargará un archivo comprimido (.zip) a su computador.
Existe tambien la posibilidad de descargar las 19 bioclimatic
variables disponibles utilizando var = ‘bio’. Por ahora
descargemos la temperatura media (tagv):
tmean <- geodata::worldclim_global(
var="tavg", # temperatura media
res=10, # resolución gruesa
path = "C:/Talleres_R/IntroR_espacial" # folder donde desea guardar los datos
)Puede ver que son 12 capas (una por cada mes, ver el atributo
nlyr en dimensions:).
## class : SpatRaster
## dimensions : 1080, 2160, 12 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## sources : wc2.1_10m_tavg_01.tif
## wc2.1_10m_tavg_02.tif
## wc2.1_10m_tavg_03.tif
## ... and 9 more source(s)
## names : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ...
## min values : -45.8840, -44.80000, -57.92575, -64.19250, -64.81150, -64.35825, ...
## max values : 34.0095, 32.82425, 32.90950, 34.19375, 36.25325, 38.35550, ...
Visualicemos dos capas (Enero y Agosto)
par(mfrow=c(1,2))
terra::plot(tmean$wc2.1_10m_tavg_01, main="Temperatura media - Enero", range=c(-50, 40), mar=c(4, 3, 4, 3), plg=list(shrink=0.9, cex=.8),
pax=list(side=1:2, cex.axis=.6))
terra::plot(tmean$wc2.1_10m_tavg_10, main="Temperatura media - Agosto", range=c(-50, 40), mar=c(4, 2, 4, 4), legend=FALSE,
pax=list(side=c(1,4), cex.axis=.6))Excelente, ya vemos algunos rásters con datos reales. Ahora tratemos de aplicar algunas funciones.
Funciones locales
Vamos a aplicar las funciones stdev()y
app() al objeto tmean que tiene 12 capas.
- Medidas de dispersion
stdev()
Veamos en donde hay mayor variación de la temperatura durante el año,
utilizando la función stdv() del paquete
terra.
tmean_sd <- stdev(x = tmean, # SpatRaster object
pop = T, # Si es verdadero (T), entonces la desviación estándard de la población es calculada
na.rm = T # Si es verdadero (T), entonces valores NA son ignorados
)¿Es lo esperado?
Tips!
Recuerde que en RStudio puede iluminar la función (ej.,
stdev) y oprimir F1 para obtener ayuda, es
decir información sobre los argumentos que se encuentran dentro de una
función.
- Calcular la temperatura media utilizando
app()
También podemos calcular la temperatura media del año, utilizando la
función app() del paquete terra.
tmean_mean <- terra::app(x = tmean, # SpatRaster
fun = mean) # función
plot(tmean_mean, main="Temperatura media anual (°C)")7.6.2. Métodos zonales y globales
Los métodos zonales y globales hacen referencia a hacer cálculos del
área total del ráster utilizando la función global() o
zonas de un ráster utilizando la función zonal() (estas
zonas se pueden definir con otro capa). Por ejemplo, la figura muestra
cómo podemos calcular la sumatoria de una zona en el ráster usando un
polígono.
Veamos un ejemplo de cómo podemos utilizar la base de datos global de la huella humana (Human footprint, en inglés) para calcular el promedio de huella humana para cada país del mundo.
La huella humana reune la mayor parte de las acciones humanas con el potencial de deteriorar la naturaleza (ver el trabajo de Venter et al 2016 para mayores detalles).
La escala es 0 = muy baja, hasta 50 = muy alta huella humana.
Note que en nuestros códigos vamos a utilizar el simbolo
%>%, conocido como forward-pipe operator en
inglés, que envia el resultado de una expresión a la otra función dentro
de un mismo objeto. Este operador%>% fue originalmente
desarrollado paquete
magrittr.
En este caso cargamos la base de datos “footprint” y después la
agregamos 10 veces utilizando la función aggregate(). Tenga
en cuenta que la función aggregate() se encuentra en varios
paquetes, entonces es buena práctica especificar que estámos utilizando
la función del paquete terra, asi:
terra:aggregate().
¿Por qué agregamos el ráster huella humana? Porque queremos la huella humana en un ráster de resolución más gruesa para acelerar el procesamiento de nuestros siguientes pasos (es decir menos tiempo de procesamiento). Más adelante veremos unos ejemplos del tiempo de procesamiento en diferentes resoluciones.
huella <- geodata::footprint(year = 2009,
path ="C:/Talleres_R/IntroR_espacial")%>%
terra::aggregate(fact=10, fun=mean)##
|---------|---------|---------|---------|
=========================================
## class : SpatRaster
## dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : wildareas-v3-2009-human-footprint
## min value : 0.00
## max value : 45.21
Ahora carguemos los países del mundo. Es una base de datos que contiene información vectorial. Más adelante dedicaremos una sección a este formato.
## class : SpatVector
## geometry : polygons
## dimensions : 231, 2 (geometries, attributes)
## extent : -180, 180, -90, 83.65625 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## names : GID_0 NAME_0
## type : <chr> <chr>
## values : ABW Aruba
## AFG Afghanistan
## AGO Angola
Ahora si, apliquemos la función zonal() del paquete
terra para calcular el promedio de la huella humana en cada
país del mundo.
huella_paises <- terra::zonal(huella,
paises,
as.raster=T, # si es F obtiene solo los datos.
fun = mean
)
plot(huella_paises)7.6.3. Seleccionado un área geográfica dentro de un ráster
Es muy común que usted tenga un ráster con una extensión global, pero
usted quiera enfocar su análisis en una región específica. Para esto
vamos a utilizar la función crop() del paquete
terra::.
Por ejemplo, seleccionemos los datos de la huella humana solo para
Colombia.
Por ahora transformemos la información vectorial a ráster, utilizando
la función rasterize() del paquete
terra::.
Note que vamos a utilizar el ráster huella como
plantilla para rasterize() los polígonos
paises; es decir obtendremos un ráster de las mismas
dimensions:, resolution:, extent
y coord. ref.:.
Para los valores del ráster utilizaremos el nombre de la columna en
el objeto vectorial que contiene el nombre de los países, es decir
NAME_0.
## class : SpatRaster
## dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## categories : NAME_0
## name : NAME_0
## min value : Afghanistan
## max value : Åland
Las categorias del objeto paises_r pueden ser
inspeccionadas con la función levels() del paquete
terra. Generalmente son dos columnas, la primera es un
numero entero que identifica los valores de la celda y la segunda es una
etiqueta que indica el nombre de la categoria. Si revisa la tabla podrá
ver que Colombia es el valor 43.
knitr::kable(head(terra::levels(paises_r)))%>%
kableExtra::kable_paper("hover", full_width = F) %>%
kableExtra::scroll_box(width = "100%", height = "200px")
|
Lo que vamos hacer es copiar el objeto paises_r y
llamarlo Colombia, luego vamos a reclasificar todos los
valores diferentes a 43 con NA (sin valores), es decir escoger
Colombia.
## class : SpatRaster
## dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## categories : NAME_0
## name : NAME_0
## min value : Colombia
## max value : Colombia
# Compare el código de arriba con este de abajo. ¿Cuales son las diferencias? ¿Podría ser de utilidad?
# Colombia <- paises_r == 43Visualicemos el mapa de Colombia. Hemos seleccionado lo
que queremos, pero el extent: sigue siendo el mismo, es
decir todo el globo entero, por eso vemos el país tan pequeño.
Podemos utilizar la función trim() para remover las
filas y columnas con valores NA (alrededor de los valores del mapa de
Colombia), y como consecuencia se cambia el extent: del
nuevo objeto.
## class : SpatRaster
## dimensions : 200, 146, 1 (nrow, ncol, nlyr)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : -79, -66.83333, -4.25, 12.41667 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## categories : NAME_0
## name : NAME_0
## min value : Colombia
## max value : Colombia
Ahora sí utilicemos la función crop() para extraer los
datos de temperatura media en el extent: del objeto
Colombia_t.
## class : SpatRaster
## dimensions : 200, 146, 1 (nrow, ncol, nlyr)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : -79, -66.83333, -4.25, 12.41667 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : wildareas-v3-2009-human-footprint
## min value : 0.00
## max value : 37.86
Visualicemos el nuevo objeto. Como podemos ver la función
crop() extrae lo valores de la huella humana con base en el
extent: del objeto Colombia_t.
Para tener solamente los valores de la huella humana dentro de los
limites de Colombia_t debemos utilizar la función
mask(). Veamos:
Ahora bien, por qué no utilizar la opción mask() desde
un principio? Hay varias razones.
- las dos rásters deben ser del mismo
extent:. Al intentar implementar la funciónmask()directamente al mapa globalhuellaobtenemos elError. Es decir el mapahuella(es global) y el mapaColombia_ttienen diferenteextent:.
## Error: [mask] extents do not match
- la segunda razón es que es más eficiente para R implementar la
función
crop(), porque selecciona filas y columnas del área definida por elextent:que queremos. Al contrario,mask()require hacer cálculos más complejos para seleccionar los pixeles que están bordeando la figura del mapa de Colombia.
Tips!
Por eso, es importante implementar primero crop() y
luego mask() para acelerar el proceso.
7.6.4. Seleccionando valores en un ráster
En ocasiones necesitamos utilizar solo unos valores determinados dentro del ráster.
Por ejemplo, queremos obtener un mapa con las temperaturas por encima
de 25 grados. Para lo cual podemos reclasificar los valores menores a 25
con NA (no values).
## class : SpatRaster
## dimensions : 1080, 2160, 1 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : mean
## min value : 25.0000
## max value : 30.9879
Pero tambien podemos necesitar un rango , por ejemplo entre 15 y 20
grados. Para eso utilizamos el operador | que indica “o”
(los valores menores a 15 “o” valores mayores a 20, serán iguales a
NA).
tmean_templado <- tmean_mean
tmean_templado[tmean_templado < 15 | tmean_templado >20] <- NA
plot(tmean_templado)Tambien podemos extraer lo valores en una pila de rásters. Por ejemplo, seleccionemos los valore mayores a 10 grados en todas las rásters.
## class : SpatRaster
## dimensions : 1080, 2160, 12 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## names : wc2.1~vg_01, wc2.1~vg_02, wc2.1~vg_03, wc2.1~vg_04, wc2.1~vg_05, wc2.1~vg_06, ...
## min values : -45.884, -44.8, -57.92575, -64.1925, -64.8115, -64.35825, ...
## max values : 10.000, 10.0, 10.00000, 10.0000, 10.0000, 10.00000, ...
7.6.5. Seleccionando valores en un ráster con otro ráster
Vamos a extraer los valores de precipitación global en areas con
temperaturas mayores a 20°C (tmean_caliente).
Primero, descarguemos la precipitación utilizando el paquete
geogata:
prec <- geodata::worldclim_global(
var="prec",
res=10,
path = "C:/Talleres_R/IntroR_espacial"
)
prec## class : SpatRaster
## dimensions : 1080, 2160, 12 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## sources : wc2.1_10m_prec_01.tif
## wc2.1_10m_prec_02.tif
## wc2.1_10m_prec_03.tif
## ... and 9 more source(s)
## names : wc2.1~ec_01, wc2.1~ec_02, wc2.1~ec_03, wc2.1~ec_04, wc2.1~ec_05, wc2.1~ec_06, ...
## min values : 0, 0, 0, 0, 0, 0, ...
## max values : 908, 793, 720, 1004, 2068, 2210, ...
Ahora calculemos la precipitación media anual:
## class : SpatRaster
## dimensions : 1080, 2160, 1 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : mean
## min value : 0.0000
## max value : 932.5833
Finalmente, utilicemos la función mask() para extraer
los valores de precipitación en areas con temperaturas mayores a
20°C.
## class : SpatRaster
## dimensions : 1080, 2160, 1 (nrow, ncol, nlyr)
## resolution : 0.1666667, 0.1666667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : mean
## min value : 0.0000
## max value : 932.5833
Visualizemos los dos rásters:
# Creemos dos columnas para los plots
par(mfrow=c(1,2))
# Temperature > 20 °C
plot(tmean_mean*0, legend=F) # Fondo gris para mostrar los continentes
t_breaks <- seq(minmax(tmean_caliente)[1], minmax(tmean_caliente)[2], 1) # escala de la temperatura
terra::plot(tmean_caliente, add=T, type="interval", plg=list(x="bottom", cex=1.2, title="°C"), breaks = t_breaks)
title("Temperatura media > 20 °C", line = -2, cex.main=1.6)
# Precipitación en áreas con temp > 20 °C
p_breaks <- seq(minmax(prec_caliente)[1], minmax(prec_caliente)[2], 150)
plot(tmean_mean*0, legend=F)
terra::plot(prec_caliente, add=T, type="interval", plg=list(x="bottom", cex=1.2, title="mm"), col=rev(topo.colors(15)), breaks=p_breaks)
title("Precipitación en áreas con temp > 20 °C", line = -2, cex.main=1.6)7.6.6. Extraer los valores de un ráster
En algún momento de sus análisis necesitará los valores de los
rásters para hacer algunos cálculos adicionales o preparar gráficas. El
paquete terra ofrece la función values() para
extraer lo valores del ráster.
Veamos un ejemplo con la temperatura media y la precipitación global.
Visualicemos una gráfica que muestre los valores pareados de temperatura media y precipitación.
Interesante ver que las precipitaciones medias anuales en áreas frías son menores que áreas más calientes.
Ahora tomemos solamente los datos de precipitación para Colombia.
Colombia_t_prec <- terra::crop(prec_mean, # precipitación media
terra::aggregate(Colombia_t, 2) # agreguemos colombia 2 veces para que sea la misma resolución de la precipitación.
)%>%
resample(aggregate(Colombia_t, 2))%>% # asegurarnos que las rásters sean iguales en extent, resolution, etc.
mask(aggregate(Colombia_t, 2)) # Seleccionar solamente el borde administrativo de Colombia
plot(Colombia_t_prec)El Pacífico siempre lluvioso!
Seleccionemos también la temperatura media para Colombia.
Colombia_t_tmean <- terra::crop(tmean_mean, aggregate(Colombia_t, 2))%>%
resample(aggregate(Colombia_t, 2))%>%
mask(aggregate(Colombia_t, 2))
plot(Colombia_t_tmean)Ahora visualicemos los valores pareados de temperatura media y precipitación para Colombia con relación al global.
Definitivamente tropical!
plot(values(tmean_mean), values(prec_mean), xlim=c(-55, 30), ylim=c(0, 1000),
ylab="",yaxt="n", xlab="",xaxt="n", pch=21, col="black")
par(new=T)
plot(values(Colombia_t_tmean), values(Colombia_t_prec), col="red", xlim=c(-55, 30), ylim=c(0, 1000), pch=22)
legend("topleft", legend=c("Colombia", "Global"), col=c("red", "black"), pch=c(22,21), cex=0.8)7.6.7. Proyectando rásters
Hasta ahora hemos utilizado coordenadas geográficas para estos análisis. Sin embargo, los análisis más locales o regionales requieren de ráster proyectadas. Es decir los datos proyectados en una superficie plana. En especial si queremos hacer cálculos de áreas para informar sobre aspectos del área de estudio (ej., área de bosques fragmentados, etc).
Hay diferentes proyecciones cartográficas, pero la que más vamos a utilizar en el análisis espacial es la que mantienen el sentido del área. Cada país y región tiene proyecciones apropiadas que preservan el area, entonces hay que hacer el ejercicio de buscar la mejor proyección para el área de estudio.
Por fortuna hay muchos sitios para buscar y spatialrefence y epsg.io.
Por ejemplo, si buscamos “Colombia” en spatialreference
nos aparecerá algo asi:
En nuestros códigos podemos utilizar cualquiera de los dos formatos ESPG (es el que aparece en el primer pantallazo de ‘spatialreference’) o Proj4.
Por ejemplo si utilizamos EPSG:3116 el equivalente Proj4
es:
+proj=tmerc +lat_0=4.5962004166 +lon_0=-74.077507916 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Veamos cómo utilizamos esto en R.
- Carguemos la cobertura vegetal (solamente la clase árboles, “trees”),
- seleccionemos los datos para Colombia y
- utilicemos la proyección
EPSG:3116para reproyectar el ráster.
La función project() de terra incluye
varios argumentos:
y =hace referencia a “coordinate reference system, (por sus siglas en inglés), en este casoEPSG:3116.method =“bilinear” es util para variables continuas y “near” para variables categóricas. Vamos a utilizar “bilinear porque es la proporción de bosque es cada celda (es decir que es una variable continua).res =hace referencia al tamaño de celda (en metros) que queremos al final. Es decir 1Km por 1Km. Celdas cuadradas, lo que facilitará y mejorará la exactitud (accuracy, en inglés) de los cálculos.
arboles <- geodata::landcover(var= 'trees',
path = "C:/Talleres_R/IntroR_espacial" )%>%
terra::crop(disagg(Colombia_t,2))%>%
terra::resample(terra::disagg(Colombia_t,2))%>%
terra::mask(disagg(Colombia_t,2))%>%
terra::project(y="EPSG:3116", method="bilinear", res = 1000)
arboles## class : SpatRaster
## dimensions : 1857, 1358, 1 (nrow, ncol, nlyr)
## resolution : 1000, 1000 (x, y)
## extent : 451351.3, 1809351, 19203.33, 1876203 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Colombia Bogota zone (EPSG:3116)
## source(s) : memory
## name : trees
## min value : 0
## max value : 1
Utilicemos la función cellSize() de terra
para calcular el área de cada celda de este ráster. Casi todas las
celdas son cercanas a 1Km2 (~1,000m x ~1,000m), pero hay algunas muy
pequeñas variaciones, en especial en el oriente de Colombia, porque esta
proyección esta centrada más en la mitad de Colombia.
Continuemos con nuestro análisis, asumiendo que las celdas son de igual tamaño. Sin embargo, tenga en cuenta reportar la exactitud en sus resultados.
7.6.8. Reclasificando rásters
La reclasificación de un ráster es de mucha utilidad en el análisis espacial. Por ejemplo si usted quiere reclasificar la cobertura de árboles en Colombia en cuatro categorias.
Para esto debemos crear una matriz que contenga el valor inicial, el valor final del ráster y que valor queremos al final. Repetir el proceso hasta terminar el rango.
Por ejemplo, definamos tres nuevas categorías de acuerdo a la proporción de árboles y reclasifiquemolas con 1, 2 y 3. En la primera linea se dice, de 0 a 0.25 reclasifiquelo a 1.
## [,1] [,2] [,3]
## [1,] 0.00 0.25 1
## [2,] 0.25 0.50 2
## [3,] 0.50 1.00 3
Ahora apliquemos esta matriz al ráster.
rc <- classify(arboles,
rclmat)
plot(rc, legend = FALSE,
col = c("dark green", "red", "yellow"), axes = FALSE)
legend("topright",
legend = c("Alta", "Media", "Baja"),
fill = c("dark green", "yellow", "red"),
border = FALSE,
bty = "n") # turn off legend borderPodríamos tener unas datos del área de cada categoría de la siguiente forma:
## layer value count
## 1 1 0 91
## 2 1 1 213491
## 3 1 2 164580
## 4 1 3 761372
7.7. Tiempo de procesamiento
El tiempo de procesamiento cuando se manejan o hacen análísis con rásters incrementa considerablemente por varias razones. Mencionemos algunas:
- La resolución espacial (tamaño de la celda o grid cell en inglés). Celdas más pequeñas generan mayor tiempo de procesamiento.
- La cobertura (tamaño de la grilla o extent en inglés). Areas más grandes requieren de mayor tiempo de procesamiento.
- La complejidad y/o tipo de la función o el análisis que se quiere hacer. Por ejemplo, reporjectar un raster puede tomar mucho tiempo.
- La combinación de todas las anteriores puede generar mucho tiempo de procesamiento o mucho consumo de memoria.
Veamos un ejemplo para visulalizar cómo se afecta el tiempo de procesamiento cuando manejamos o hacemos cálculos con rásters.
Resolución espacial
Creemos 3 grillas de la misma extensión pero diferente tamaño de celda:
pixel <- c(1, 10, 100, 1000, 5000, 7500, 10000) # Tamaño de celda
pixel_size <- list()
for(i in 1:length(pixel)){
pixel_size[[i]] <- rast(ncol=pixel[i], nrow=pixel[i])
values(pixel_size[[i]]) <- 1:ncell(pixel_size[[i]])
}Ahora apliquemos una función sencilla a cada ráster y calculemos el tiempo que tarda cada una en procesar.
tiempo_proceso <- list()
tiempo <- list()
for(j in 1:length(pixel_size)){
tiempo_proceso[[j]] <- system.time(sqrt(pixel_size[[j]]))[3]
}##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
## pixel tiempo
## 1 1 0.00
## 2 10 0.00
## 3 100 0.00
## 4 1000 0.01
## 5 5000 0.22
## 6 7500 2.17
## 7 10000 3.92
Visualicemos en una gráfica para que pueda observar la tendencia.
Por ahora tengamos en cuenta que el tiempo de procesamiento es una limitante cuando manejamos rásters. Las soluciones a este problema van desde la optimización de sus códigos hasta la utlización de supercomputadores para acelerar los procesos.
7.8. Tamaño de los rásters creados por R
Cuand se carga un ráster o se aplica un función a un raster, R trata de utilizar la memoria RAM de su computador. Sin embargo, si los valores del el objeto (es decir un ráster) no pude ser guardado en la memoria RAM, entonces R empieza a crear archivos temporales que generalmente son bastante grandes.
Los archivos temporales se crean cuando no se asigan la ruta y nombre del archivo dentro de la función o en funciones en las cuales no existe esta opción de asignar la ruta y el nombre al archivo final.
Antes de ver unos ejemplos, recordemos que al principio habíamos configurado un folder para guardar los archivos temporales.
## [1] "C:/Talleres_R/IntroR_espacial/Temp_para_borrar/r_tmp_2024-07-21_195156.559877_28580_30265.grd"
Ahora si veamos unos ejemplos. Carguemos un raster (elevación)
utilizando el paquete geodata::. La función
elevation_global() le permite indicar la ruta
if(!dir.exists("C:/Talleres_R/IntroR_espacial/Input_data")){
dir.create("C:/Talleres_R/IntroR_espacial/Input_data")
}else{
print("Folder Input_data ya existe")
}## [1] "Folder Input_data ya existe"
## class : SpatRaster
## dimensions : 4320, 8640, 1 (nrow, ncol, nlyr)
## resolution : 0.04166667, 0.04166667 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : wc2.1_2.5m_elev.tif
## name : wc2.1_2.5m_elev
## min value : -415
## max value : 7412
Ahora intentemos desagregar este ráster, es decir
## [1] 373248
## [1] 37324800
## --- none ---
##
|---------|---------|---------|---------|
===========
## [1] 149299200
## --- none ---
## [1] "C:/Users/jzulo/AppData/Local/Temp/RtmpMVZkU4/spat_UCzrRLqRud378rP_28580.tif"